clc; clear; close all

% Importing the temperature data. The process is same for all 
% the dataset. Before importing, save the source file as a copy
% and delete all the header leaving just the temperature data
% at first row. 
% This program will give you heatflux output in "Heat flux.txt" 
% file format. 
% While allocating the temperature_T1 to T4, the thermocouple
% placements might have changed. Here T1 is the largest temperature of four
% that T1,T2,T3 and T4 and will decrease in that order.
% So, the proper arrangement for the temperature profile is required.
% Figure 2 will provide the temperature fitting and it should have a
% decreasing slope which is obtained by ensuring T1 being largest 
% temperature value and T4 being the smallest one. 

temperature_input = importdata('Enter the temperature file (.lvm) here'); 
temperature_t = temperature_input(:,1);
temperature_T1 = temperature_input(:,5);
temperature_T2 = temperature_input(:,2);
temperature_T3 = temperature_input(:,4);
temperature_T4 = temperature_input(:,3);


figure
subplot(2,1,1)
plot(temperature_t,temperature_T1,'Color',[0 0 0],'LineWidth',1);
hold on
plot(temperature_t,temperature_T2,'Color',[0.4 0.4 0.4],'LineWidth',1);
hold on
plot(temperature_t,temperature_T3,'Color',[0.6 0.6 0.6],'LineWidth',1);
hold on
plot(temperature_t,temperature_T4,'Color',[0.7 0.7 0.7],'LineWidth',1);
hold on

ylim([80 300])
xlabel('Time, \itt\rm (s)')
ylabel('Temperature, \itT\rm (^\circC)')
set(gca,'FontName','Arial')
set(gca,'FontSize',16)


k = 392; % Thermal conductivity of Cu, W/mK.
tc_location = [0 2.54 5.08 7.62] * 1e-3;
temp_profile = [temperature_T1 temperature_T2 temperature_T3 temperature_T4];
n = 4;
slope_denominator = n * sum(tc_location .^ 2) - (sum(tc_location)).^2;
slope = (n * sum(tc_location .* temp_profile,2) - sum(tc_location) .* sum(temp_profile,2)) ./ slope_denominator;
q = - k * slope / 1e4; % heat flux, W/cm^2;

subplot(2,1,2);
plot(temperature_t,q,'LineWidth',1);
xlabel('Time, \itt\rm (s)')
ylabel('Heat flux, \itq\rm (W/cm^2)')
set(gca,'FontName','Arial')
set(gca,'FontSize',16)

[Q, id_max] = max(q);

figure
scatter(tc_location,temp_profile(id_max,:),'sb');
P = polyfit(tc_location,temp_profile(id_max,:),1);
temp_fitted = P(1) .* tc_location + P(2);
hold on
plot(tc_location,temp_fitted,'r--');
xlabel('Thermocouple location, \itz\rm (m)')
ylabel('Temperature, \itT\rm (^\circC)')
set(gca,'FontName','Arial')
set(gca,'FontSize',16)

Ts = P(1) .* 13.1825 * 1e-3  + P(2);



disp(Ts);
disp(Q);

writematrix(q,'Heat flux.txt');
